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Abstract. A new slender-body theory for viscous flow, based on the concepts of dimensional 
reduction and hyperviscous regularization, is presented. The geometry of fiat, elongated, or 
point-like rigid bodies immersed in a viscous fluid is approximated by lower-dimensional objects, 
and a hyperviscous term is added to the flow equation. The hyperviscosity is given by the 
product of the ordinary viscosity with the square of a length that is shown to play the role of 
effective thickness of any lower-dimensional object. Explicit solutions of some simple problems 
illustrate how the present method is able to represent with good approximation both the velocity 
field and the drag forces generated by rigid motions of the immersed bodies, in analogy with 
classical slender-body theories. This approach has the potential to open the way to more 
effective computational techniques, since the complexity of the geometry can be significantly 
reduced. This, however, is achieved at the expense of involving higher-order derivatives of the 
velocity field. Importantly, both the dimensional reduction and the hyperviscous regularization, 
combined with suitable numerical schemes, can be used also in situations where inertia is not 
negligible. 



1. Introduction 



Composite systems where microscopic bodies move in a viscous fluid are ubiquitous in both 
biological and technological contexts, and the flows involved are often characterized by very small 
Reynolds numbers. In an early seminal contribution, Stokes ( 1851 ) computed the flow past a 



translating rigid sphere in the low- Reynolds-number limit, determining the corresponding drag 
force and showing how Newtonian liquids can override D'Alembert's paradox. This result was 



surrounding flow was provided by Kim & Karrila ( 2005 ) 



extended to ellipsoidal particle shapes by Oberbeck ( 1876 ), and to the case of rotating and shearing 



surrounding flows by Edwardes (1892) and Jeffery (1922), respectively. The solution for a general 



An explicit solution of the same problem for a rigid body of general shape immersed in a low- 
Reynolds-number flow appears to be out of reach, but the broad spectrum of possible applications 
have favoured the development of methods to approximate the velocity field generated by the 
body, or at least the drag force and torque exerted on it. For a rigid body whose length is large 



compared with its breadth, starting with an idea of Burgers ( 1938 ) , various authors, including Tuck 



(1964), Tillett (19701, Batchelor (1970), and Cox (1970), have contributed to the development of 



what is known as slender-body theory. In essence, that theory provides approximate expressions 
for the relevant quantities simulating the action of the immersed body on the surrounding fluid 
by means of force and torque distributions concentrated on suitable lines. The flow generated 
by those distributions is divergent precisely on the lines where forces concentrate, hence the use 
of the name "singularity method" to describe this technique. The same methodology has been 



exploited also for particles which are not slender since the formative c ontribution of |Oseen ( 19271 
and a clear account of it with applications to Stokes flows is given by Chwang & Wu ( 1975 1 . 



Despite the possibility of approximating with good accuracy many quantities of physical in- 
terest, it is difficult to apply classical slender-body theory to particles of arbitrary shape and to 
situations featuring the presence of many particles, which are clearly important for applications. 
For this reason, modifications of the theory have been proposed by de Mestre (1973), Johnson 

and |Barta (2011), and also somewhat different theories have been 



(1980), Barta & Liron 



developed by Lighthill 



(1975 19761 and Keller & Rubinow (1976). In a recent paper, Cortez & 



Nicholas (20121 presented a new theory in which the concentrated force densities are replaced by 
forces which are localized in a small three-dimensional region lying within the slender body. Our 
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approach can be viewed as complementary to theirs, in that we consider not only concentrated 
force distributions, but also concentrated bodies, regularizing in turn the flow equation. 

We propose a new slender-body theory which can be applied to rigid bodies whose slenderness 
can affect one dimension (flat bodies), two dimensions (elongated bodies), or even three dimensions 
(point particles). In so doing, we are motivated by the understanding that the hyperviscous 
regularization of the Stokes equation described in £j2] makes it is possible to obtain a solution for 
the flow past a translating point particle such that: (i) the drag exerted by the particle on the 
fluid is equal to the total viscous traction exerted on a Newtonian fluid by a spherical particle; 
(ii) the solution of the regularized problem is a good approximation of the classical solution 
for a slender sphere. In <|3j we derive the Green's function for the regularized Stokes operator, 
which is exploited, in Sj4j to compute the flow generated by a point particle, and to evaluate the 
hydrodynamic interaction between point particles. The problem of the flow generated by a rigid 
rod is considered in Ej5j where some limitations of the singularity method, associated with the 
rigidity constraint, are considered. In the final discussion, we point out possible advantages in 
numerical simulations based on our slender-body theory, and we outline further developments. 



2. Hyperviscous regularization 



The classical Stokes equation for incompressible Newtonian fluids, that is the low-Reynolds- 
number linearisation of the Navier-Stokes equation, reads 



du 



= —Vp + pAu + pb , 



where p is the pressure field, u is the divergence-free velocity field, p > is the constant and 
homogeneous mass density, p > is the dynamic viscosity and pb is a volumetric force density. 
For a steady flow, the Stokes equation reduces to 

Vp — pAu = pb . 

It is well-known that the regularity of solutions to the three-dimensional Navier-Stokes equation 
is still an open problem, and various modifications of that equation with better regularity theories 



have been analysed. Among those, the hyperviscous regularization (see, for example, |Lions (19691, 
Chap. I, Remarque 6.11) consists of adding a term proportional to AAu to the equation. For this 
modified equation, the existence and uniqueness of regular solutions (that is of solutions which are 
continuous on the flow domain for every instant in a finite time interval) have been established. 
The low-Reynolds-number linearisation of that equation is, for a steady flow, 

Vp - pAu + £AAm = pb , 

where the additional parameter £ > is called the hyperviscosity. 

Despite the mathematical appeal of the hyperviscous regularization, it is generally hard to 
assign a relevant physical meaning to £, other than that of a higher-order dissipation coefficient. 



In a series of papers, Fried & Gurtin (20061, Musesti (2009), Giusteri et al. (20111, and Giusteri 



(2012) introduce and analyse different contributions to £ associated with dissipation functionals. 



Here, we wish to assign to £ a geometric, rather than dynamical, meaning, which allows for a 
clearer interpretation of the new parameter in terms of standard physical quantities. To this end, 
we introduce a new length-scale L > 0, and set £ = pL 2 , so that the hyperviscous steady flow 
equation becomes 

Vp- pA(u- L 2 Au) = pb. (1) 



2.1. Dimensional reduction and hyperviscosity. In the slender-body theory presented in this 
paper, L is viewed as the effective thickness of the lower-dimensional objects which approximate 
three-dimensional slender bodies. We now clarify how the hyperviscous regularization discussed 
above can be combined with dimensional reduction to approximate the flow past slender bodies, 
and the role of L in this procedure. 
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Recall the Stokes problem of determining the steady flow v past a rigid sphere S of radius L 
centred at the origin, 

' Vp v - fiAv = in R 3 \ S, 
divv = in E 3 \ S, 

(2) 

v = on <9E, 
v = U at infinity, 

where U £ K 3 denotes the constant far-field velocity. The hyperviscous Stokes problem for the 
steady flow u past a point particle placed at the origin has the form 

Vp u - fiA(u - L 2 Au) = dS(x) in M 3 , 

divw = inM 3 , (3) 

u = at the origin, 

where 8 is the three-dimensional Dirac distribution and d is the drag force. 

The solution to the Stokes equation in the entirety of M 3 is divergent at the origin; the solution 
of (|3| (with L > 0) is bounded and continuous in 0, making explicit the reason why the hypervis- 
cous equation is considered a regularization of the classical one. It is important to notice that, in 
the limit L — > 0, the problem ^ becomes ill-posed. 

A far-field condition akin to that entering |2| is not imposed in ([3|, but the latter problem 
admits a unique solution, and the value of u at infinity is determined by the force d. In {4.1 we 
solve explicitly problem ([3]), showing that, far enough from the sphere of radius L, u is a good 
approximation of the solution v of problem Moreover, u takes at the origin the same value 
taken by v on the spherical surface <9£, and it turns out that, to obtain also u = U at infinity, we 
must take the drag force exerted by the point particle as 

d= -[i \ (Vv + Vw f ) n = 6TTfj,LU , 

where n denotes the unit outward normal to <9£ and Vu' denotes the transpose of Vd. Namely, 
d must be equal to the traction exerted on the fluid by the sphere in problem ([2|. It is precisely 
this last fact which suggests the interpretation of L as the effective thickness of the point particle. 

A point particle represents a body which is slender in every spatial direction, but we are 
interested also in describing the motion of flat as well as elongated bodies. Actually, there are 
two different classes of problems worthy of mention: resistance problems and mobility problems. 
In a resistance problem, the velocities of the immersed objects are assigned and the goal is to 
determine the drag force required to sustain the motion. In a mobility problem, the forces acting 
on the immersed objects are prescribed and the goal is to determine the resulting velocity field. 

The solution of resistance and mobility problems in the case of a point particle provides a 
linear relation between the assigned velocity and the drag force, which is the only force acting 
on the point particle. The situation is remarkably different when dealing with flat or elongated 
rigid bodies. Whereas in resistance problems the prescription of the velocity of an immersed body 
automatically includes the rigidity constraint, the same constraint can be imposed in mobility 
problems only by assigning as force density the sum of the drag force and the reactive forces that 
are generated within an immersed body as a consequence of the assumption that it is rigid. We 
comment more on this in f|5j 

3. Green's function for Stokes flow 

The basic tool used to construct solutions to the Stokes problem is the Stokeslet, that is the 
Green's function for the Stokes operator in M 3 . Its classical expression is given by 

h ■ x 

h (h-x)x 

s(x) = S(x)h = — + HW , 

y J V ' 8tt^i\x\ 87r/i|a;| 3 ' 
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where h £ M 3 is any fixed vector and S is the Oseen tensor, with Cartesian components 



Sij(x) 



87r^i|a;| 8TTfi\x\ 3 ' 

where Sij is Kronecker's symbol. 

To compute the regularized Stokeslet, we first need a Green's function g solution of the fourth- 
order elliptic equation 

fiL 2 AAg - fxAg = S(x-x'). 
Using the Fourier transform, we easily obtain 

d£ 3 (k). 



g(x - x') = 



(27T) 



fi\k\ 2 {L 2 \k\ 2 + 1) 

We choose a base for the momentum space in such a way that x — x' is along the ^-direction, set 
R = \x — x'\, switch to polar coordinates (fc, 0, (j>), and use the calculus of residues to obtain 



g(x - x') 



2tt 



-oo pi ^ikRcosl 



(2ir) 3 fiL 2 J J_j k 2 + l/L 2 



d(cos 0) dk 



sin kR 



(2ir) 2 vL 2 Rj k(k 2 + l/L 



(2tt) 2 hL 2 R 



Im 



+oo 



2 ) 

ikR 



dk 



k(k 2 + l/L 2 ) 



dk 



(2ir) 2 f iL 2 R 



ttL - ttL e L 



Hence, the desired Green's function is 
g(x - x') = 



1 — exp — 



\x — x 

IT" 



4tt/i\x — x' 

Notice that, in the limit L — >• 0, Q reduces to the fundamental solution 

1 



gi(x - x 1 ) 



Ann\x — x'\ 



(4) 
(5) 



for the Laplace equation. 



3.1. Regularized Stokeslet. We now proceed to construct the regularized Stokeslet, that is a 
pressure field and a velocity field £ satisfying 

divC = 0, (6) 
Vp c + AC = hS(x) , (7) 

with /i£l 3 and 

A = [iL 2 AA - fiA . 

Let 4> satisfy A(f> = 5(x — x'). Then, since A commutes with V, a solution for |6])-([7]) is given by 

PC = -Ad , 
C = hcf> + W. 

The scalar field $ entering this solution is chosen to satisfy the constraint ([6| and turns out to 
have the explicit form 

where g\ is as defined in ^ and the convolution /i * fi of two scalar fields /i and f% is defined by 

(fl*f 2 )(x) := J h{x-x')f 2 (x')dx' . 

Now, exploiting the properties of the convolution and the operator A, and using the Green's 
function g given by equation Q, we find that 

* A{h ■ Vg)) = -{gi * div(A(gh))) = - div(. gi * A{gh)) = h • V gi 
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and 



51 * (ft ■ Vg) = 



1 



i(h-fc)5i5e lfc - x d£ 3 (fc) 



(2tt) 3 
ft ■ x 



1 



i(h • k)e 



ih'X 



(2tt) 3 J ^|fc| 4 (i 2 |fe| 2 + l) 

+ °° rl i cos 6*e ifc l x l cos6 

d(cos 0) dk 



dC 3 {k) 



4ir 2 pL 2 \x\ J Q J_ 1 k{k 2 + l/L 2 ) 
sin(fc|a;|r) 



ft • x 



+ 00 



An 2 fiL 2 \x\ J_ t J k{k 2 + l/L 



dk dr 



ft ■ x 

87r^i|a;| 
ft ■ x 



rl (1-e- l |t| ) dr 



1 2L i.i 



87r/i|a; 

Hence, our regularized Stokeslet is given by 



2L- i _ M , 
e t — 1 



R 2 v £ 



ft, • a; 

47r|a;| 3 ' 



C(») 



ft 



87r/i|a; 



1 - 2e 



2L 



x 



2L 2 



x 



(ft • x)x 



1 + 2e" 



6L 



6L 2 / _ 



87T^i|a;| 3 

Notice that the regularization does not alter the pressure field. Moreover, the velocity field £ 
is well defined for any x G IR 3 , as is the fundamental solution <?, at variance with the classical 
expression s, which is singular at the origin. 

We also define the tensor field Z, the regularized Oseen tensor, such that, using Cartesian 
components, 



Z l3 {x) 



87r^|a:| 



1 - 2e 



2L 



2L- ( 



(.-¥-!) 



+ 



X 1 X n 



l + 2e 



6L 



M 6L 2 ( _w 

L + \xT {< ~ 



8TTfl\x\ 3 

whereby it follows that C( x ) — Z(a;)ft. The Stokeslet allows us to obtain an integral representation 
for the solution of (JlJ, with vanishing condition at infinity, in the form of a convolution: 



u(x) :=p I Z(x - x')b(x') dC 3 (x') . 
Jr 3 



(8) 



Due to the regularity of the regularized Stokes tensor Z at the origin and the properties of 
convolution, u as determined by pi) is bounded and continuous on all of space, even if the force 
density pb is localized on a set of zero Lebesgue measure — in other words, even if 6 is a measure 
concentrated on a lower-dimensional subset of M. 3 , (JsJ) determines a well-behaved velocity field. 
This makes it possible to model interactions between the fluid and lower-dimensional objects by 
assigning force distribution concentrated on them, and thereby to turn the approximate singularity 
solutions of the classical theory into exact solutions of an approximated model. Notice that 
although ([8| readily offers a solution for mobility problems, the force density b must represent the 
interactions taking place in all of space, not only in the fluid region; hence, knowledge of both the 
drag force applied on the rigid body and the reactive forces generated within the body is required. 

4. Flows generated by slender spheres 

4.1. A single point particle. The mobility problem for a point particle subject to a drag force 
density dS(x) (viewed in a reference frame in which the point particle is at rest) can be readily 
solved on applying formula ([8]), and gives the velocity field 

u(x) = Z(x)d - Z(0)d = Z(x)d - (QirpL^d . (9) 
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The resistance problem requires a bit more of work. Recall that the disturbance field generated 
by a sphere rigidly translating with constant velocity U in a Newtonian liquid is given by 



U , U (U-x)x ,3(U-x)x 

u (x) = -a— - b— - - a " ' h 6 — 

v ' \x\ \x\ 3 \x\ 3 x 5 



(10) 



where the constants a and b can be determined by imposing the no-slip condition u s = U on 
the surface of any sphere of finite radius. However, since u s is divergent at the origin, we cannot 
impose the no-slip condition on a sphere with vanishing radius, namely a point particle. 

For a hyperviscous liquid, we seek the divergence-free disturbance field u produced by a trans- 
lating sphere. It satisfies the equation 

Vp - fiA(u - L 2 Au) = . (11) 

Since, given a velocity field u s defined on all of space, it is always possible to find a pressure field 
p s such that Vp s — fiAu s = 0, any solution of 

u- L 2 Au = u s (12) 



provides a solution for (11), with pressure p s . It is readily checked that u 9 = u s + L 2 Au s is a 



particular solution of ( 12 ). We write the general solution of the associated homogeneous equation, 
using spherical coordinates (r, #,</>), as 



+00 n 



«° = E E e n > l hW(ir/L)Yl(6,<(>) 



n—0 l— — r 



where c™' ( are vectors of coefficients, hn are spherical Hankel functions of the first kind, and 
are spherical harmonics. Any solution of ( 11 ) is then of the form u = u 9 + u°, and the boundary 
conditions determine the coefficients. 

We now wish to show that it is possible to obtain a disturbance field that is bounded at the 
origin and satisfies the no-slip condition on a sphere with vanishing radius. We assume without loss 
of generality that U = U^e^, and we write the Cartesian components of u in spherical coordinates 
(the angular dependence of being understood) , obtaining: 



ui(r,M) = l/3 



36 + 6aL 2 



u 2 (r,0,0) = ill 3 



3b + 6aL 2 



-f-OO 



(r+^ + EE^^^)^ 



n=0 l=—n 



J7 3 2aV167r 
3~r 



U 3 {3b + 6aL 2 



16tt 



+oo n 



^ + EE C 3' i ^ 1) ^/ £ ) y n 



n—0 l— — n 



Thanks to the orthogonality of the spherical harmonics and the properties of the spherical Hankel 
functions, imposing iti(0, 9, </>) = 7/2(0, 9, cf>) — 0, we find, for any value of a, that 6 = and that 
the other non-vanishing coefficients are given by 

2.1 2.-1 . 2.-1 . 2.1 2U3CL 2ir 



L V 15 

The third condition, namely it 3 (0, #,</>) = U3, implies that C3' = U^^/Att, c 2 ' = U^^/ 'An /5, and 
a = —3L/4, while all the remaining c^' 1 vanish. 
Hence, the disturbance field is given by 



ui{x) 



3L U3X3X1 



l + 2e~ T 



6L 



6L 2 



(.-*-.)■ 
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u 2 (x) = 



3L U 3 x 3 x 2 



1 + 2e~ 



61 



x 



6L 2 



t - l 



w 3 (a;) 



1 - 2e 



M 2L 
\xf 



lasl 2 



3LU 3 x 



- 1 



1 + 2e~ 



6L 



6i 2 



(.-*-!) 



We thus find that the velocity field, in a reference frame in which the point particle is at rest, is 

u(x) = 6w(iLZ(x)U- U. (13) 

Comparing ^ and (13), we conclude that the drag force needed to sustain the latter motion is 
given by the celebrated Stokes relation d — 6irfj,LU. Finally, when \x\ > L, the difference between 
the velocity field given by (13) and the classical solution of ([2]), obtained by (10) with a = — 3L/4 
and b = — L 3 /4, is small and decays as \x\~ 3 in the far field. 

4.2. Two point particles. For the resistance problem related to the motion of two point particles 



translating with constant velocity, it is tempting to exploit the correspondence developed in S4.1 



and solve instead the mobility problem, using ^ and the superposition principle. That procedure 
works only when the two point particles move with the same velocity, because the flow produced 
by the motion of the particles must be steady. 

Let us consider the disturbance field w produced by the two particles translating with velocity 
U. We apply the drag force density 

pb(x) = di5{x - y{) + d 2 5(x - y 2 ) 

to the system, where y 1 and y 2 are the positions of the particles in a co-moving frame. Now the 
solution of the mobility problem reads 

w(x) = Z(x - y x )di + Z(x - y 2 )d 2 , 

and, to find the solution for the resistance problem, we must be able to choose d in such a way 
that w(y 1 ) — w(y 1 ) = U. We have the conditions 

Z(0)di + Z( Vl - y 2 )d 2 = U = Z{0)d 2 + Z(y x - y 2 )d x , 

which readily imply d\ = d 2 =: d, and 

d=[Z(0) + Z(y 1 -y 2 )}- 1 U. 

Without loss of generality, we may take one of the coordinate axes to be parallel to y x — y 2 , and 
we can easily see that d is parallel to U only when the latter is either parallel or orthogonal to 
V\ ~ Vi- I n an y other case, we find that a component of the force which is not in the direction of 
the motion is needed. That additional component of force is required to cancel the hydrodynamic 
interaction which tends to move the points away from the desired trajectory, since the leading 
point (the one in front with respect to the direction of U) pulls and is pushed by the other point 
with a force directed along y x — y 2 . 

4.3. Hydrodynamic interaction between point particles. Even though the above techniques 
cannot provide the full solution to the unsteady situation in which two point particles move with 



different velocities, it is possible to utilize a procedure, known as method of reflections (Kim & 



Karrila 2005), to approximate the instantaneous flow field and the instantaneous hydrodynamic 



interaction between the particles. So far, we have considered particles moving in a fluid which is 
at rest at infinity, but it is easy to include the effect of more complicated conditions in the far- 
field, as described by a field w^. Indeed, the disturbance field produced by a particle uniformly 
translating with velocity U in an ambient flow is 

u(x) = 6npLZ(x - y)[U - w 00 (y)} , 

where y is the position of the particle in a co-moving frame. 



G. G. GIUSTERI AND E. FRIED 



The method of reflections for two particles (labelled 1 and 2) consists of an iterative scheme 
which, at every step k, takes as ambient field for a particle tOoo plus the field generated by the 
other particle in the step k — 1. Now we have i,j = 1,2, i ^ j, and 

uf\x) = 6tt[iLZ(x - y^iUi - w^y^} , 
=Z{x- yi )F?\ 

where U\ and U2 are the instantaneous velocities of the two particles, and 
Ff ] = en^LiU, - Woo ( Vl ) - Z( Vi - yj )Ff- l) ) . 

Then we wish to express F j 2fe ' and F[ 2k+1 ^ in terms of 

F<f } =6TTtiL(U 1 -w 00 (y 1 )) 

and 

F ( 2 ] = 6irfiL(U 2 ~ w 00 (y 2 )) . 
After some manipulations, we obtain 

fe— 1 

p (2k) = Q 2k p(0) + (q^ F (0) + Q2n+1 p(0)\ (w) 

n=0 

and 

k— 1 

F (2k+l) = F (0) + Q 2fe+1 F (0) + (q 2 «+2 F (0) + Q2n+1 F (0)\ ^ ^ 

n=0 

where Q :— —Gtt^lL Z{y 1 — y 2 ). 

Notice that, with every application of Q, an additional factor of L/\y 1 — y 2 \ appears. It therefore 
becomes easy to keep track of the order of the approximation in terms of that small parameter: 
the A:-th approximation F^ is exact up to the order (L/\y 1 — y 2 \) k ~ 1 - The same is not true if we 
extend the procedure to the case of N point particles, because of the appearance of corrections to 
lower-order terms even after many iterations. Indeed, let N identical point particles be placed at 
distinct points y i moving with velocity Ui, i = 1, . . . , N. Taking i, j = 1, . . . , N, we have 

«i ( x ) = 6tt[iLZ(x - yi )[Ui - WooiVj,)] , 
F 4 (0) =&KtxL{U i - Woo {y i )), 
u<*>=Z(*- Vi )lf\ 



where we have set 



F\ k) = 67r M i(c/ i - - E z (y« - yj) F . 



(fc-i) 

i 



Considering now indexes io, i\, ■ . ■ = 1, . . . , N , we have 

fe 



10 10 Z-^i 



e n f 



(0) 



i\ ,...,i n \m— 1 



(16) 



where Qkj — Qjk '•— —SjkGTtfJ,LZ(yj — y fc ), and the product in the above formula is the standard 
multiplication of matrices. Though somewhat complicated, expressions (14 1, (15), and (16) can 
be used to approximate the hydrodynamic interactions in dilute suspensions of slender spheres, 
that is suspensions in which the average separation between particles is much greater than L. 
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5. Approximating a rigid rod 



We now consider the case of a rigid straight rod, represented by a line segment of length 2a, 
lying on the xi-axis in the interval [—a, a], when seen in a co-moving frame. We are interested in 
the steady disturbance field produced by a constant rigid motion of the rod, represented by the 
rigid velocity field 

w(x) = Qx + v , 

where Q is the skew-symmetric spin (namely, a constant and uniform spin tensor), and v is a 
constant translational velocity satisfying Q.v — 0; the latter compatibility condition forces v to be 
parallel to the axis of rotation, and is required to ensure a steady flow. More explicitly, the spin 
matrix takes the form 

o -n 3 n 2 ' 
n= n 3 o -fix 

-J? 2 fix 

where fli is the spin about ej. We work in a co- moving frame and, as a consequence, the flow at 
infinity matches —w. In this setting, we may assume that the steady Stokes flow has the form 
u — w, where u represents the disturbance flow, and is such that u = w on the rod and u = at 
infinity. 

We are again facing a resistance problem; however, we first attempt to solve a mobility problem, 
introducing a force density concentrated on the rod A := [—a, a] x {(0,0)} in the form 

pb(x) = f{xx)S A (x) = f(xx)x[-a,a](xx)S(x 2 )S{x 3 ) , 

where X[-a,a]{ x i) is the characteristic function of [—a, a], and 5a represents a uniform probability 
measure supported on A. Using the convolution (|S|) , we have 



2 hWx) 



1 - 2e"T 



2L 



2L 2 



-e L 



E 



M " M 2 

a Qi ~ SuXj)(xj - Sxjx'^fjjx'x) 
8tt M M3 



e L 



dx'x 



1 + 2e~T 



6L 



6L 2 
M 2 



(e-* - l) 



dx'x, 



where we have introduced yj (xx — x'-^) 2 + x\ + x 2 =: M. We evaluate the previous expression at 
x 2 = x 3 = 0, giving 



Ui (aii, 0,0) = 



fijxx - y) 

8nfi\y\ 
+ Sx 



„ |y| 2L 

1 - 2e"T _ e - 
1 2/1 



Xl+a fi(xx~y) 
Xl - a 87rju|y| 



2L 2 / 










L 


-o; 








ivi 6L 


Ivl 


6L 




L H e 


L 






|y| 




+ M 





dy. 



Since we are interested in the solution for |rpx| a : to impose the matching with Wi on A, we may 
split the integral in two parts, uf and u~ , accordingly evaluating |y|: 



u+{xx) 



u t (xx) 



Xl+a h(xx-y) 

8irny 
+ Sx 

-Mxx-y) 
8^2/ 



2L 2 2L^_ 



1 - 2e^r 

y z y 
Xl+a h(xx-y) 



2L 2 



dy 



o 

l-2e* 
o 



5xi 



87r/xy 

2L 2 2L 

H 5- H e 

-fi(xx - y) 



1 + 2e" 



2L 2 

y 2 



-e £ 



dy 



6L^ 

,,2 



6L 2 



-e t 



„ v 6L y 6L 2 

1 + 2eT e i =- 

2/ 2T 



6L 2 



dy. 
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In agreement with the classical slender-body theory, proposed, for instance, by Batchelor (1970), 
we make the following Ansatz on the concentrated force field: 



fi(r) = SiT + c.i with Si, Cj € R. , for i = 1, 2, 3. 
In the sequel, we will need the following indefinite integrals, for ael and n G 



(17) 



r 1 ™ _1 

/ ^ ay dy = -^ v H 
j y fe=i 



(n-l-fc)!^" 1 a' 
(ra-1)! y n - k + Jn- 1)! 



rEi(ay) ; 



here, as usual, Ei denotes the exponential integral function, defined as 

Ei(y) = / —dr. 



We also recall the series expansion of Ei around y — 0, which has the form 

y k 

k\k 7 



Ei(y) = 7 + log |y| + ^ V 



fc=i 

where 7 is the Euler-Mascheroni constant. 

We next evaluate and u~ for a force of the form (17) and find that they are given by 

— I — / \ Ci 1 SiX± 

uj{x x ) = 



87T/X 



1 (jxt+a + §u j% 



87T/i 



Jo 1+a + SuJo 



and 



with 



u, (xi) 



-(Ci + SiXi) / o 



87T/i 



X 



7-^1+0 



logy-Ei(-2//i) + -e * + ^ (e * - l) 



3L 



I? 

>r 

3L 2 



x\-\-a 



logy-Ei(-y/L) e l-^-fe - 1 

y r 



(-*-') 



S x\—a 

y° = 



2L 2 

y + 2Le~% H ( e~i - 1 



(.-*-!) 



or -* 6L " 1 - - , 
y — 2Le l e & — 1 

y 







log I ?/ 1 - Ei(y/L) ei 

y 

3L 



y 2 

3i 2 



ei -1 



log|y|-Ei(y/L) + — er 

y y 



(.* - 1) 



:ci —a 




/v xi — a 



y - 2Le T 
y + 2Le* - 



2L 2 

y 

6L 2 

y 



e* -1 



(.*-!) 



a;i — a 




a;i — a 



where evaluation at stands for the limit, as necessary. Now, we have 
Ui(xi,0,0) = uf(xi) + u~(xi) 



Cj + SjXi fTXl + a x f Xl +a , ; o r -f)0 \ 

— — — (i. +ou-L - y xi _ a - ouy xl - a ) 



8717/ 



(18) 
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We now take Sj = 0, as classical procedures would suggest that, to sustain a rigid translation, 
we should compute the velocity in the middle of the rod, equate it to U , and thereby arrive at 
the following approximation for the constant force densities, valid to any order in e = L/2a: 

2~K[LV\ 



h{x'x) = c 2 
h(x'i) = c 3 



log^" 1 ) + 7 - 1/4 + £ 2 

AlTfjLV 2 

log(e- 1 )+ 7 + 2-e 2 ' 

47r/i«3 
log(e- 1 ) + 7 + 2 - e 1 ' 



(19) 
(20) 
(21) 



These expressions agree with results obtained by Batchelor (1970), but, unfortunately, the distur 



bance field given by ( 18 ) is not constant on A, and hence u ^ U except at the midpoint of the 
rod! 

Indeed, such a solution is almost constant in a neighbourhood of the origin, and it would still 
be possible to impose a genuinely constant velocity on the boundary of a slender body of non- 
vanishing breadth, when A lies well within it. This justifies the classical approximation for slender 
bodies, but, since our object is one-dimensional, we really need to obtain a velocity field which is 
exactly constant on the whole of A. As already mentioned, the previous discrepancy arises from 
the fact that we must include in the force density / the reactive force which guarantees the rigidity 
of the rod. An analogous result is obtained if we set Ci = 0, to study a rigid rotation. 

We could find a hint about the form of such forces by actually solving the resistance problem 
for a rigidly translating rod, applying the procedure of §4.1| Indeed, with A being the degenerate 
fundamental ellipse of a family of prolate spheroids, we may view it also as a degenerate prolate 
spheroid. The solution of this problem using prolate spheroidal coordinates and external spher- 
oidal wave functions is provided in Appendix [XJ but that solution provides an expression for the 
disturbance field which is not readily expressed as a convolution. To determine the reactive forces, 
we therefore rely on a different method. 

For a fixed N € N, we approximate the rod by a sequence of equidistant points with first 



components {x \} 
component 



chosen in [—a, a] so that x\ = — a and 



2 N +1 

pb 2 (x) = S(x 2 )S(x 3 ) J2 
k—l 



= a. We then find the second 



of the force density required to make all the points move with equal velocity U = U 2 e 2 . (The 



same procedure produces analogous results for a velocity U = U\e\. 
1, 2, . . . , 2 + 1, can be found by solving the Toeplitz system 



The coefficients /g , k = 



a 


a i 


"2 


a3 


' ■ a 2 « 


Ox 


Of) 


ax 


a 2 '•■ 




«2 


ax 


a 


ax '■■ 




«3 


(to 


ax 


a 





Q,2 N 

with a, = Z 22 (ia/2 N - 1 ,0,0), for i 
two point particles. 

We solved the system for N = 



a i 



a i 
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u 2 




fi 




u 2 




f 2 N +l 
L J 2 




_u 2 _ 



0, 1, . . . , 2 , which generalizes the technique used in ( 4.2 



for 



1,2,..., 13, fixing a — 1cm, L/2a = 0.005, U 2 — 1 cm/s, 



and /i = 4 Pas. Plots of the maximum and minimum coefficients, f 2 and f 2 



against N are 



provided in Figure[T] to show that the result approaches an asymptote for N sufficiently large. The 
profile of the force density distributed on the rod, calculated for a, U 2 , and as above and different 
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Figure 1. Force coefficients fl (left) and / 2 (right) for different values of TV 
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Position 
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Figure 2. Force density distributed on A for different values of L/2a (log-scale). 



values of L/2a by using 8193 points (N = 13), is shown in Figure [2] It is much larger (by two orders 
of magnitude) towards the ends of the rod, and rather flat in the middle, as expected, since the 
motion of the central points receives the maximum support by the hydrodynamic force generated 
by the motion of the other points, whereas the endpoints receive the least contribution from such 
interactions. We see that, in accordance with the interpretation of L as effective thickness, the 
magnitude of the force needed to sustain the motion increases with L. 



6. Summary and perspectives 

We have presented a new slender-body theory, which, in contrast to previously developed alter- 
natives, might be viewed as a modelling scheme, based on the concepts of dimensional reduction 
and hyperviscous regularization. We propose to approximate the geometry of flat, elongated, or 
point-like rigid bodies immersed in a viscous fluid with lower-dimensional objects, while adding a 
hyperviscous term to the flow equation. The hyperviscosity is given by the product of the ordinary 
viscosity with the square of a length that replaces the characteristic size of the body along the 
dimensions that are shrunken to zero in the slender-body limit. 
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We explicitly solved some simple problems to show that the proposed method is able to represent 
with good approximation both the velocity held and the drag forces generated by rigid motions 
of the immersed bodies, in analogy with classical slender-body theories. The regularized Stokeslet 
and the regularized Oseen tensor have been computed. These objects are basic to the solution 
of mobility problems, in which the concentrated force density acting on the system is prescribed. 
We explained the limitations of this kind of solution technique, which stems from the need to 
determine the reactive force arising from the rigidity constraint. That force can be computed 
starting from the solution of resistance problems, in which the velocity of the immersed objects is 
prescribed. 

The construction of a numerical scheme for the simulation of slender bodies in viscous fluids is a 
necessary step for any realistic application, and will be the object of a forthcoming paper. Within 
the present approach, the computational strategy stands to benefit from the understanding that 
the unsteady movement of the slender bodies is no longer represented by a time-varying domain, 
but rather by a time-dependent constraint on the function space of admissible velocity fields. The 
complexity of the geometry is reduced at the expense of involving higher-order derivatives of the 
velocity field in the flow equation, and a time-dependent constraint on it. We nevertheless believe 
that this approach has the potential to enable more effective computational techniques, especially 
when dealing with suspensions featuring large numbers of particles. 

Another important aspect of our approach is that both the dimensional reduction and the 
hyperviscous regularization, combined with suitable numerical schemes, can be used also in situa- 



tions where inertia is not negligible. Indeed, as Giusteri et al. (2010) show, the flow generated by 



one-dimensional rigid bodies, moving in a fluid governed by a generalization of the Navier-Stokes 
equation with hyperviscosity, can be uniquely determined. A generalization to deformable bodies 
also seems feasible. 

The first author would like to express his gratitude to the Department of Mechanical Engineering 
of the University of Washington for the hospitality offered during the development of this research. 



Appendix A. Uniform flow past a straight rod 

To compute the low-Reynolds-number flow past a uniformly translating rigid rod, we follow 
the method of {4.1 Specifically, we start from the solution u ps , given by Chwang & Wu (1975), 
for the disturbance field generated by a prolate spheroid with axis along e\ and foci at x\ = ±a, 
uniformly translating with velocity U = U\e.\ + U 2 e 2 . In view of the symmetry of the spheroid, 
there is no loss of generality in taking U3 = 0. A particular solution of ( 12 ) is u* = u ps + L 2 Au ps . 
Setting 



«3 1 



Ri := yj{xi + a) 2 +r 2 . 



D n := 



R2 ~ {x\ - a) 
R\ - (xi + a) 



D, := 



1 



1 



m m 1 

D 5 := 

we can express its components as 



Di := 



D 4 := 



R> 



1 

xi - R 2 



R 2 := y / (x 1 - a) 2 + r 2 
x\ + a x\ — a 



Do := 



Ri 

x\ — R\ 



Ri 



Rj-R 2 (x 1 -a) Rf-Ri(x 1 +a) 



R\ 



1 

#2(2:1 



- 1{a x 4 
+ a\D 2 



h)logD -2P 1 x 1 D 1 

X2 



R\ - +a) ' 

h 2L 2 a 1 x 1 D 3 + {fi 2 x 2 - 2/3 1 x 1 )D 4 



a 2 x 2 Di + j3 2 -j (Ri - R 2 



+ 



x 2 xi(xi + a) (2L 2 a 2 (3 2 
~RT~ ~ R2 



2L 2 a 2 ^D 1 
x 2 xi(xi - a) ( /3 2 

Ri 



2L 2 a 2 
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u* 2 (x) = (j3 2 - a 2 ) log D - (aix 2 + 2/3 1 x 2 )D 1 + 2L 2 a 1 x 2 D 3 



+ 1 2L2a2 



xi + a ( 2L 2 a 2 
p 2 tt 2 



R 2 J r 2 V Ri 
4L 2 \ 2L 2 a 2 x 2 ( x\ + a x\ — a 



R\ 



Rl 



+ 



foxl 



1 



2fi 2 



+ (/3 2 x 2 - 2 / 3 1 x 1 x 2 ) J D 5 , 



u^(x) = - (aix 3 + 2/3 1 x 3 ) J D 1 + 2L 2 ai x 3 D 3 



02X2X3^2 



+ 



/3 2 X 2 X 3 



(.ti - a) 



1 2i?A / 1 

2- J - (si + a) 



2L 2 a 2 X2X 3 



(xi + a) 



i?i r 

1 2 



2ii 2 



i?2 



(/3 2 x 2 x 3 - 2/3iXix 3 )D 5 , 



where ai, a 2 , Pi, and /3 2 are constants to be determined. 

Now we introduce prolate spheroidal coordinates (u, v, </>), defined by 



x\ = a cosh u cos v . 



x 2 = a sinh u sin v cos < 



x 3 = a sinh u sin v sin < 



with u G [0, +00), v G [0, 7r], and G [0, 2tt). The surfaces given by u — c, with c > 0, are confocal 
prolate spheroids with foci at x\ = ±a, and u = describes the degenerate spheroid A. Using 
these coordinates, we have 

r = a sinh m sin v , R\ = a(coshw + cosw) , R 2 = a(coshw — cosu) , 



i? 2 — (xi — a) = a(cosh?i + 1)(1 — cosw) , i?i — (xi + a) = a(coshw — 1)(1 — cos v) , 
and we can rewrite the components of u* as 



*/ 1 \ r,/ o \ i ( coshu + 1 

u\{u,vA) = -2 ai+A log - 

\ cosh u — 1 



cosh w cos v cosh u sin w 
4 \osh^-cos 2 , +2ai ^ 



cosh m — cos 2 v 



4L 2 ai cosh u cos 2 1>(3 cosh 2 u + cos 2 v) cos </> sinh u sin w cos w 

^ 2 ; T~2 o — 2a 2 -g 

a (cosh u — cos z vy cosh u — cos z i> 

cos cos w 4L 2 a 2 cos cos i; 

+ 2p 2 ^-, ; h 



+ 



+ 



sinh u sin v ' a 2 sinh u sin v (cosh 2 u - cos 2 t>) 
cos 4> cosh m cos u(cosh u cos v + 1) / 2aL 2 a 2 a(3 2 
sinh u sin v \ i?f i?2 

cos 4> cosh m cos w(coshu cos v — 1) / a/3 2 2aL 2 a 2 



sinh u sin w \_ i?i i?| 

+ a(/3 2 sinh u sin w cos — 2fi Y cosh u cos v)D 4 , 
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u%(u,v,4>) = ((3 2 - a 2 )lo 
+ 



' cosh u + 1 



- (2c*i +4/?i) 



\ cosh u — 1 

4i 2 a 1 cos sinh it sin 7; cos u(3 cosh 2 u + cos 2 u) 



cos </> sinh it sin v cos v 
cosh 2 u — cos 2 v 



a 2 (cosh^ u — cos 2 v) 3 

cosh u cos v — 1 / f5 2 Ri 2L 2 a 2 \ cosh u cos v + 1 ( fi 2 R 2 2L 2 a 2 



sinh 2 u sin 2 - 



ai? 2 



sinh 2 u sin 2 u 



+2a 2 



cos 2 6 cosh 71 sin 2 v 



cosh 2 u — cos 2 v 



4L 2 



+/3 2 cos 2 



(xi - a) 



a 2 sinh 2 77 sin 2 v 
2R X 



2L 2 a 2 cos 



aRi 

xi + a x\ 



Rl 



1 



-Ri a 2 sinh 2 it sin 2 v 



(x 1 + a) 



2R 2 



1 



R 2 a? sinh 2 u sin 2 - 



+(/?20 2 cos 2 4> sinh 2 77 sin 2 w — 2j3\a 2 cos cosh u cos 1/ sinh 71 sin v)D$ , 
sin (A sinh u sin D cos v 



u* 3 (u,v,4>) = -(2a 1 + 4(3 1 )- 



cosh 77 — cos 2 v 



AL oi\ sin </> sinh u sin 7; cos v(3 cosh 77 + cos v 
a 2 (cosh 2 77 — cos 2 v) 3 

1 2i?i 



2a 2 



cos 4> sin </> cosh 77 sin 7; 
cosh 2 it — cos 2 v 



+{3 2 cos</>sin0 (xi — a) 
+2i 2 a2 cos <j> sin 0(xi + a 
— 2i 2 a 2 cos tfi sin </>(xi — a 



-Ri a 2 sinh 2 u sin 2 w 
1 2 



(xi + a) 



2i? 2 



1 2 sinh 2 u sin 2 tj 



-R? i?i<z 2 sinh 2 u sin 2 7; 
1 + 2 
-R2 i?2a 2 sinh 2 u sin 2 w 



+(/32a 2 cos (j> sin </> sinh 2 it sin 2 u — 2/?ia 2 sin </> cosh 77 cos u sinh 77 sin u)-Ds , 

Notice that, in the limit u — > 0, u* is divergent. Since we wish to impose the velocity of the fluid 
precisely on the set A defined by 77 = 0, we must add a suitable solution of the homogeneous 
equation associated with (12 1, thereby cancelling out the divergent terms. 



We introduce external prolate spheroidal wave functions (using the notation of Erdelyi et al. 



( 1955 )) defined for n G N and I = 0,...,n: 

W± l (u,v,(f>) = S 1 ^ (cosh 77, K 2 a 2 /4)Ps^ (cos v, K 2 a 2 /4)e ±ll,p . 

Those are solutions of AW + k 2 W — 0, regular at infinity but divergent on A. In our problem 
K = i/L, and we are interested in the cases n = 0, 1, 2, 3. We have 



W °(u,w» = S , ° (3) (coshu,-a 2 /4? 2 )Ps!](cosi;,-a 2 /4/ 2 



2L 



■ cosh 1 



'K. 



2r>0 2p>Q 

W?(u,v,(/>) = S , 1 0(3) (coshu,-a 2 /4? 2 )Ps?(cosw,-a 2 /4/ 2 ) 



! cosh 1 



P 2 ° (cosw) , 



»2 E E 



(-i) p < r < P ^ 



2L 



■ cosh 



2r>-l 2p>-l 

Wf^u,^) = 5 1 1(3) (coshit,-a 2 /4Z 2 )PsJ(cosi>,-a 2 /4Z 2 )e ±# 



i cosh 
L 



'J R 2 p+i(cosw) , 
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W${u,v,<f>) = ^ (3) (coshu,-a 2 /4? 2 )Ps°(cosw,-a 2 /4/ 2 ) 



2r>-2 2p>-2 v 

W^{u,v,<j>) = S , 2 1(3) (coshu,-a 2 /4? 2 )Ps2(cosw,-a 2 /4Z 2 )e ±,;0 



e ±i<p I 1 



'2 

1 



cosh 



2 



E S i-^S^ (^) *«<~" 



X 

2r>-l 2p>-l 



W|= 2 (u,«,0) = (cosh u,-a 2 /4p)Ps2 (cos u,-a 2 /4Z 2 )e d 



-.■ ' i ', u-. 

cosh u , 



2r>0 2p>0 v 7 

°( 3 )^c,l,.. 2//U2-vp o 0/„„ „„, _„2/„,2^ 



W^ 3 u (m,w» = S , 3 W (coshu,-a 2 /4? 2 )Ps^(cosw,-o 2 /4/ 2 ) 



2r>-3 2p>-3 v 



W^{u,v^) = S , 3 1(3) (coshu,-a 2 /4/ 2 )PsJ(cosw,-a 2 /4/ 2 )e ±i<# ' 



cosh u , 



X 

2r>-2 2p>-2 



VK 3 ±2 (u,w» = 53 (3) (coshu,-a 2 /4Z 2 )Ps 2 (cosw,-a 2 /4/ 2 )e ±i2 * 

= e ±i ^ (l VW'x 

V cosh u J 

x E E(-l)X>i r /^^«»A- 2 ,. + ,(^)p|„ t 3(co 8 »), 

2r>-12p>-l v 7 

where sjj and ajj m are coefficients dependent on the parameter k 2 o 2 /4 = — a 2 /4Z 2 , K u is a modihed 
Bessel function of the second kind, and is an associated Legendre function of the first kind. 

The divergence of the logarithmic terms in u* and u 2 can be cancelled adding multiples of Wq . 
If U 2 = 0L2 = P2 — 0, which corresponds to a uniform translation of the rod along its axis, the 
dependence on <j) and the radial decay of u* suggest that a solution u" of (12 1 can be written in 
the form 

u\ = u{ + Cl W$ + diW? + k{W% + e{W% , 

u\ = u* 2 + c 2 {Wl + W^ 1 ) + d 2 (W 2 1 + W 2 X ) + k 2 {Wl + W3 1 ) , 

4 = ul + ic 3 (W£ - W{ 1 ) + id 3 (W% - W 2 V ) + ik 3 {Wl - W^ 1 ) . 

The no-slip condition on A requires that u(0,v,(f>) = —U, which fixes the value of a\, once the 
boundedness of u" is ensured by a suitable choice of the the coefficients Cj, di, hi, and e\ as 
functions of a±. 
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0, corresponding to a translation of the rod 



of ( 12 I of the form 



On the other hand, in the case Ui = cti = /3j 
with direction in the plane orthogonal to e±, we may seek a solution u 

uj; =ul+ Cl {Wl + Wi 1 ) + dx{Wl + W 2 ~ x ) + k x {Wl 
u\ =u* 2 + c 2 W* + d 2 (W 2 + W 2 ~ 2 + A 2 W° + B 2 W°) 

+ fc 2 (W| + W 3 2 + A 3 W° + B 3 W°) , 
4 = u* 3 + id 3 (W 2 - W 2 2 ) + ik 3 (W 2 - W 3 2 ) . 

As above, the no-slip condition on A fixes the value of a 2 , once the boundedness of u 1 - is ensured 
by a suitable choice of the coefficients Ai, Bi, Ci, di, and ki. Finally, the disturbance field u for a 
uniform translation with generic velocity U = U\e\ + U 2 e 2 is given by u = it" + u 1 - . 

Notice that, since we can impose the corresponding boundary conditions on the surface de- 
scribed by u = c, for any value c > 0, the previous expressions provide the general solution for the 
flow past any translating spheroid of the family with foci at x\ = ±a. 
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